setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
file=read.csv("20201120/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",head=T)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

test2=file[file$DC.FDR.sig==12,]
test2=test2[test2$pattern.nodupl.num.DC==1,]

site.22=test2
site.4=read.csv("20201216/4.Coordinates.site.csv",header = T)

#site.4=tidyr::separate(site.4,unitID,into=c("chr","pos","ref","alt"),sep=":")
#site.4$unitID=paste(site.4$chr,site.4$pos,site.4$ref,site.4$alt,sep=":")
site.4$id=paste(site.4$chr,site.4$pos,sep=":")
site.4=filea[filea$unitID %in% site.4$unitID,]

#AJHG
ajhg=read.table("E:/0 公共数据库差异情况/db_for_5hmc/2704_ASMs_hg18ToHg19_AJHG.bed",head=F,sep="\t")
ajhg$id=paste(ajhg$V1,ajhg$V2,sep=":")

a=intersect(as.character(site.22$id),as.character(ajhg$id))
b=intersect(as.character(site.4$id),as.character(ajhg$id))

#imprinted loci
library(ChIPpeakAnno)
imprint=read.table("E:/0 公共数据库差异情况/db_for_5hmc/imprint.hg19verion.bed",head=F,sep="\t")
names(imprint)=c("chr","start","end")
bed.imprint =toGRanges(imprint,format="BED",header=TRUE)

bed.site.22=data.frame(chr=site.22$Chr,start=site.22$Start-1,end=site.22$Start)
bed.site.22=toGRanges(bed.site.22,format="BED",header=TRUE)

bed.site.4=data.frame(chr=site.4$Chr,start=as.numeric(site.4$Start)-1,end=as.numeric(site.4$Start))
bed.site.4=toGRanges(bed.site.4,format="BED",header=TRUE)

ol1=findOverlapsOfPeaks(bed.imprint,bed.site.22)
t=as.data.frame(ol1$overlappingPeaks)

ol2=findOverlapsOfPeaks(bed.site.4,bed.imprint)
t2=as.data.frame(ol2$overlappingPeaks)

write.csv(t,"./20210127.117K.vs.22/22loci.overlap.imprinted.loci.csv",quote=F,row.names = F)


#ASOC
library(openxlsx)
asoc.iPSC=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",sheet=6,startRow = 2)
asoc.iPSC$id=paste(asoc.iPSC$Chr.x,asoc.iPSC$Start.x,sep=":")
asoc.NPC=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",sheet=7,startRow = 2)
asoc.NPC$id=paste(asoc.NPC$Chr.x,asoc.NPC$Start.x,sep=":")
asoc.iN.Glut=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",sheet=8,startRow = 2)
asoc.iN.Glut$id=paste(asoc.iN.Glut$Chr.x,asoc.iN.Glut$Start.x,sep=":")
asoc.iN.GA=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",sheet=9,startRow = 2)
asoc.iN.GA$id=paste(asoc.iN.GA$Chr.x,asoc.iN.GA$Start.x,sep=":")
asoc.iN.DN=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",sheet=10,startRow = 2)
asoc.iN.DN$id=paste(asoc.iN.DN$Chr.x,asoc.iN.DN$Start.x,sep=":")

asoc.iPSC = data.frame(id=asoc.iPSC$id,group=rep("iPSC",length(asoc.iPSC$id)))
asoc.NPC = data.frame(id=asoc.NPC$id,group=rep("NPC",length(asoc.NPC$id)))
asoc.iN.Glut = data.frame(id=asoc.iN.Glut$id,group=rep("iN.Glut",length(asoc.iN.Glut$id)))
asoc.iN.GA = data.frame(id=asoc.iN.GA$id,group=rep("iN.GA",length(asoc.iN.GA$id)))
asoc.iN.DN = data.frame(id=asoc.iN.DN$id,group=rep("iN.DN",length(asoc.iN.DN$id)))

asoc=rbind(asoc.iPSC,asoc.NPC,asoc.iN.Glut,asoc.iN.GA,asoc.iN.DN)
asoc.id=unique(c(as.character(asoc.iPSC$id),as.character(asoc.NPC$id),as.character(asoc.iN.Glut$id),
                 as.character(asoc.iN.GA$id),as.character(asoc.iN.DN$id)))
length(intersect(as.character(site.22$id),as.character(asoc.id)))
length(intersect(as.character(site.4$id),as.character(asoc.id)))

#TWAS

twas.BD=read.table("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/BD_Ruderfer2018.dat",head=T,sep="\t")
twas.BD=twas.BD[!stringr::str_detect(twas.BD$TWAS.P,pattern = "NA"),]
twas.BD$TWAS.P=as.numeric(twas.BD$TWAS.P)
twas.BD$TWAS.FDR=p.adjust(twas.BD$TWAS.P,method = "BH")
twas.BD=twas.BD[twas.BD$TWAS.FDR<0.1,]

length(intersect(unique(unlist(strsplit(site.22$Gene.refGene,";"))),twas.BD$ID))#0 OVERLAP
length(intersect(unique(unlist(strsplit(site.4$Gene.refGene,";"))),twas.BD$ID)) #1 OVERLAP

TWAS.SCZ.BD=read.table("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/BDSCZ_Ruderfer2018.dat",head=T,sep="\t")
TWAS.SCZ.BD=TWAS.SCZ.BD[!stringr::str_detect(TWAS.SCZ.BD$TWAS.P,pattern = "NA"),]
TWAS.SCZ.BD$TWAS.P=as.numeric(TWAS.SCZ.BD$TWAS.P)
TWAS.SCZ.BD$TWAS.FDR=p.adjust(TWAS.SCZ.BD$TWAS.P,method = "BH")
TWAS.SCZ.BD=TWAS.SCZ.BD[TWAS.SCZ.BD$TWAS.FDR<0.1,]

length(intersect(unique(unlist(strsplit(site.22$Gene.refGene,";"))),TWAS.SCZ.BD$ID))#0 OVERLAP
length(intersect(unique(unlist(strsplit(site.4$Gene.refGene,";"))),TWAS.SCZ.BD$ID)) #3 OVERLAP

TWAS.ASD=read.table("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/AlzheimersProxyMetaIGAP_Marioni2018.dat",head=T,sep="\t")
TWAS.ASD=TWAS.ASD[!stringr::str_detect(TWAS.ASD$TWAS.P,pattern = "NA"),]
TWAS.ASD$TWAS.P=as.numeric(TWAS.ASD$TWAS.P)
TWAS.ASD$TWAS.FDR=p.adjust(TWAS.ASD$TWAS.P,method = "BH")
TWAS.ASD=TWAS.ASD[TWAS.ASD$TWAS.FDR<0.1,]

length(intersect(unique(unlist(strsplit(site.22$Gene.refGene,";"))),TWAS.ASD$ID))#0 OVERLAP
length(intersect(unique(unlist(strsplit(site.4$Gene.refGene,";"))),TWAS.ASD$ID)) #3 OVERLAP

TWAS.Depr=read.table("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/Depression_Nagel2018.dat",head=T,sep="\t")
TWAS.Depr=TWAS.Depr[!stringr::str_detect(TWAS.Depr$TWAS.P,pattern = "NA"),]
TWAS.Depr$TWAS.P=as.numeric(TWAS.Depr$TWAS.P)
TWAS.Depr$TWAS.FDR=p.adjust(TWAS.Depr$TWAS.P,method = "BH")
TWAS.Depr=TWAS.Depr[TWAS.Depr$TWAS.FDR<0.1,]

length(intersect(unique(unlist(strsplit(site.22$Gene.refGene,";"))),TWAS.Depr$ID))#0 OVERLAP
length(intersect(unique(unlist(strsplit(site.4$Gene.refGene,";"))),TWAS.Depr$ID)) #3 OVERLAP
